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ABSTRACT 

We investigate the stability of gravitational accretion of an ideal gas onto a compact object moving 
through a uniform medium at Mach 3. Previous three-dimensional simulations have shown that such 
accretion is not stable, and that strong rotational 'disk-like' flows are generated and accreted on 
short time scales. We re-address this problem using overset spherical grids that provide a factor of 
seven improvement in spatial resolution over previous simulations. With our higher spatial resolution 
we found these 3D accretion flows remained remarkably axisymmetric. We examined two cases of 
accretion with different sized accretors. The larger accretor produced very steady flow, with the 
mass accretion rate varying by less than 0.02% over 30 flow times. The smaller accretor exhibited an 
axisymmetric breathing mode that modulated the mass accretion rate by a constant 20%. Nonetheless, 
the flow remained highly axisymmetric with only negligible accretion of angular momentum in both 
cases. 

Subject headings: accretion — hydrodynamics — shock waves 



1. INTRODUCTION 

The gravitational accretion of an ideal gas onto a 
compact object has been applied across a wide range 
of astrophysical topics, including the growth of super- 
massive black holes an d associated feedback in AGNs 
(|Booth fc Schavdl2009t ). accretion onto young star clus- 
ters as an expl anation for the forma tion of ultraluminous 
x-ray sources (|Naiman et al.l 12009ft an d accretion onto 
gas g iant planets around evolved stars (jVillaver fc Livid 
2009). A common application is that of wind-fed X- 
ray binaries, beginni n g wit h the original suggestion by 
iDavidson fc Ostrikerl |l973| ) that the X -ray luminosity 
of Cyg X-3 is driven by the gravitational capture of ma- 
terial in the stellar wind of the companion star. More 
recently, wind-fed accretion onto neutron stars has been 
invoked to explain the behav ior of supergiant fast x-ray 
transients (jDucci et al.H2009l) . 

While the basic theory for gravitational accretion onto 
a moving star was d e scribe d more than 70 years ago 
bv iHovle fc Lvttletonl (|1939f ). from a theoretical stand- 
point we know surprisingly little about its stability and 
temporal behavior. The simplified case of planar ac- 
cretion in two dimensions is known to be unstable to 
the growth of small pe rturbations from a steady state 
(jBlondin fc Pope 1 12009ft . leading to accretion through a 
quasi-Keplerian accretion disk. Whether similar behav- 
ior exists in three dimensions remains an open question. 
Previous 3D simulations have observed rotational flow 
around the surface of the accretor, but were limited by 
the computational resources necessary for small accre- 
tor radii and sufficiently high resolutions. In this paper, 
we describe the results of high-resolution hydrodynamic 
simulations that probe the question: is three-dimensional 
Hoyle-Lyttleton accretion stable, or is the axisymmetry 
of the flow broken through rotational flow or violent flip- 
flop behavior? 

A thorough review of Hoyle-Lyttleton accretion (HLA) 
theory and subsequent numerical studies is given by 
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lEdgarl (|2004D . We provide only a brief synopsis here. 
Hoyle and Lyttleton derived a formula for the accretion 
radius of a point mass M moving at a speed Voo through 
a uniform medium of density pea'- 



2GM 



(1) 



Gas approaching the star with an impact parameter 
less than R a will collide on an accretion line behind the 
star and will be left with insufficient kinetic energy to 
escape the gravitat i onal p otential of the star. From this, 
IHovle fc Lvttletonl (|1939ft posited that the mass accre- 
tion rate is given by the mass flux through a circle of 
radius R a far upstream: 



(2) 

There have been many numerical simulations of HLA 
in 2D, be ginning with t he axisymmetric steady -state so- 
lutions of iHuntl (11971ft . iMatsuda et~aH (|1987ft were the 
first to observe what t hey called the 'flip-flop' instability, 
which was also seen bv lFrvxell fc Taaml Jl988). Flip-flop 
behavior is characterized by an accretion wake that oscil- 
lates between alternately rotating disk states with a brief 
transitory phase in which the disk is flushed into the ac- 
cretor. This non-axisymmetric instability has been sub- 
sequently seen in numerous hydrod ynamic simulations of 
two-dimensional planar accretion (Matsuda et al. 1991; 



Zarinclli et al. 1995; Bcncnsohn et al. 1997; Shim a et al.l 
19981: IPoeorelov et all 12000ft . IBlondin fc Pope I (|2009ft 
showed that the flip-flop instability was a true oversta- 
bility rather than a consequence of numerical effects, and 
that the growth rate increases as the radius of the accre- 
tor is decreased. 

The relevance of instabilities in two-dimensional HLA 
to the temporal behavior of accreting X-ray sources has 
been debated since the discovery of the flip-flop insta- 
bility. Given that accretion-powered X-ray sources in 
Be-type binary systems are expected to accrete from a 
planar disk expelled by the companion Be star, attention 
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has focused on the relationship of the flip-flop inst ability 
to these Be-type X-ray pulsars. iTaam et"aLl ([1988D inter- 
preted the recurrent flares in EXO 2030+375 in terms of 
the ep isodic accretion r e sultin g from the flip-flop insta- 
bility. iBenensohn et al.l (|1997| ) argued that the variable 
accretion torques associated with the flip-flop instability 
may explain the erratic spin behavior in these systems. 

It remains to be seen whether this disk accretion mode 
exists in 3D. The best simulations to date of 3D HLA 
were a series of studies using neste d cart esian grids 
(jRuffert k Arnetti [l99l IRuffertl UM% [1991 . The ac- 
cretion in these 3D studies was found to be unstable, 
although not as violent as the flip-flop behavior found 
in 2D plana r flow. Nonethele s s, som e of the 3D models 
reported in iRuffert k Arnettl (|1994 hereafter RA) did 
exhibit a quasi-periodic swinging of the accretion shock 
leading to brief epochs of strong rotational flow near the 
accrctor. This behavior was particularly evident for the 
models with the smallest accretors. The models with an 
accretor radius of R s = 0.05 R a (S6 and S8) exhibited 
rotational flow near the surface of the accretor, but the 
flow velocities were sub-Keplerian by more than a factor 
of two and the rotational flow pattern lasted for no more 
than 10 orbits. 

However, it is important to note that these simula- 
tions were of relatively low spatial resolution and limited 
to relatively large radii for the accreting star. Although 
these 3D runs used nested cartesian grids to improve the 
resolution near the accreting star, each nested grid was 
only 32 zones on a side. Their best resolution amounted 
to 6.4 zones across the radius of the accretor. This is very 
low compared to that used in many of the 2D simulations 
exhibiting flip-flop behavior (Blondin k Pope 2009 used 
a resolution 6 times finer). Moreover, these simulations 
were limited to an accretor radius of 0.05 R a (and one 
short run at 0.01 R a ) due to limited computational re- 
sources. This limitation is important because the growth 
rate of the flip-flop instability in 2D i ncreases dramati- 
cally with decreasing accretor radius (|Blondin k Popel 
2009), with long phases of quasi-disk accretion only oc- 
curring for much smaller accretors (0.0037 R a )- 

In this paper we describe two (R s /R a = 0.05 and 0.01) 
high-resolution simulations of three-dimensional HLA. 
The parameters of these models we re ch osen to mimic 
two of the simulations described by IRAI in which they 
found time-dependent, non-axisymmetric behavior. 

2. COMPUTATIONAL METHOD 

The simulations presented here use the VH-1 hy- 
dr qdynamics code, w hich is the same code described 
in iBlondin k Pope I ((2009). VH-1 uses a Lagrange- 
remap formulation of the p iecewise parabolic method 
(|Colella k Woodward! [1981 to solve the Euler equa- 
tions for the inviscid flow of a compressible ideal gas. 
The piecewise parabolic method is an explicit Godunov 
method that offers high-order accuracy (3rd-order in 
space, 2nd-order in time). 

The extension to 3D is enabled by the use of a Yin- 
Yang overset grid as described by iKagevama k Satol 
(2004) . Given that one needs to follow the accretion flow 
over orders of magnitude in distance from the central ob- 
ject, a spherical coordinate system is the natural choice 
for accretion studies. Unfortunately in 3D this choice 
introduces a coordinate singularity associated with the 



polar axis that creates severe problems with the Courant 
limit on the time step in explicit methods. It also results 
in unwanted artifacts associated with advection near the 
axis. The Yin- Yang approach solves these problems by 
using two spherical polar grids aligned by 90 degrees with 
respect to each other. The polar regions of one grid are 
replaced by the equatorial region of the other, with over- 
lapping regions used to implement boundary conditions 
in the angular directions. 

We use Yin and Yang grids with dimensions of 
252x72x216 zones, which when combined result in a grid 
circumference of 288 zones. We use a non-uniform ra- 
dial grid, keeping Ar/r roughly constant throughout the 
grid. This grid enables us to simulate the downstream 
behavior of the accretion shock extending to 10i? s while 
still offering a resolution near the accrctor comparable to 
that of the highest-resolution 2D simulations published 
to date. The number of radial zones was chosen to pro- 
duce a value of Ar/r ss 0.021, comparable to the angu- 
lar spacing of 7r/144. This ensures approximately cubic 
zones throughout the grid. For the R s — 0.01 R a run we 
extended the radial grid to 324 zones to keep the same 
value of Ar/r. 

We orient our grid so that the incoming gas flows to- 
ward the origin from the +x direction, such that the 
projected face of the Yin- Yang grid is laterally symmet- 
ric (Figure [J). This places the accretion wake almost 
entirely in the Yang grid, and minimizes the impact of 
the grid seams on the subsonic post-shock flow. During 
preliminary runs, we observed the formation of small, 
low-entropy streamers that originated in the low-density 
stagnation region of the post-shock flow and propagated 
upstream toward the accretor. These streamers would 
occasionally wrap partly around the accretor, disturbing 
the flow slightly as they were accreted and causing minor 
motion of the wake. The origin of these streamers ap- 
pears to lie in the seams between the Yin and Yang grids, 
as reorienting the flow so that the seams are no longer 
present in the accretion wake causes them to vanish. 




Figure 1. A schematic diagram of the grid orientation, as seen 
looking down the +z axis. The incoming flow begins on the Yin 
grid (purple) and flows toward the +x axis. The bow shock (green) 
is almost entirely contained within the Yang grid (blue). The grid 
extends to the outer edges of the highlighted plane. A portion of 
the grid has been made semi-opaque in the middle to illustrate the 
Yin-yang grid geometry. (A color figure is available in the online 
journal.) 



The inner boundary of our simulation represents the 
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surface of the accretor. In order to minimize the im- 
pact of this inner boundary on the accretion flow, we 
imposed an absorbing boundary condition by setting the 
density and pressure to very small values inside of this 
boundary. This is simila r to the inner boundary cond i- 
tion used bvlHuntl (11971ft. iRAl . iBenensohn et all (|1997ft . 
and lBlondin fc Pope I (|2009V To allow gas to flow unim- 



the steady bow shock. 



peded off the grid downstream of the accretor, we imple- 
ment a zero-gradient outer boundary downstream, but 
with the added condition that gas may not flow inward 
onto the grid once it has escaped. 

In order to monitor the accretion of mass and angu- 
lar momentum, we record the total mass and angular 
momentum crossing the inner boundary over each time 
step. By averaging these values periodically over small 
intervals of time, we are able to obtain a measure of the 
mass and angular momentum accretion rates. 

To initialize the 3D grid, our approach is to begin with 
2D axisymmetric simulations on a spherical grid. These 
2D simulations are evolved for several flow times to allow 
them to settle into a steady state (if such a state exists), 
after which the evolved 2D solution is mapped onto the 
3D Yin- Yang grid. The 3D grid uses the exact same 
radial grid as in the 2D simulation, avoiding the need 
for interpolation in the radial direction. Moreover, the 
2D and 3D codes share the exact same FORTRAN code 
for hydrodynamic evolution. This procedure is meant 
to minimize any transients generated by the initial con- 
ditions on the 3D grid. Starting from a steady state 
allows us to investigate the origin of any instability by 
watching small perturbations grow over time. It also al- 
lows us to avoid the relatively long initial evolution that 
previous researchers had to compute while waiting for 
large-amplitude transients to die away. 

Simulation parameters are chosen to match models S8 
(S6 was the same model but lower resolution) and T10 of 
IRAl . with an upstream Mach number of three, a ratio of 
specific heats of 5/3, and accretor radii of R s = 0.05 R a 
and 0.01 R a . Assuming a strictly uniform plane-parallel 
flow at the outer boundary has the adverse effect of cou- 
pling the behavior of the flow at the accretor boundary to 
the outer radius of the grid because the outer boundary 
is not truly an infinite distance away. To more closely 
match the true HLA problem, we calculate the velocity, 
density, and pressure (assuming isentropic flow) at the 
outer edge of the grid by solving for the ballistic trajec- 
tory from infinity ( Bisn ovatvi-Kogan et al.lll979ft . 
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Figure 2. The mass accretion rate as a function of time for the 
model with R s /R a = 0.05. 



This steady flow retains the original axisymmetry to a 
high degree. The maximum specific angular momentum 
of the accreted gas never exceeds 0.05% of the specific 
angular momentum of a Keplerian orbit at the radius of 
the accretor. This axisymmetry is illustrated in Figure 
[3] which shows the gas pressure in a 2D plane normal to 
the incoming flow but downstream of the accretor. 




3. RESULTS 
3.1. Large Accretor: R s /R a = 0.05 

With the l arger accretor (R s /R a = 0.05, corresponding 
to S8 in IRAl ) . we observed a remarkably stable flow over 
60 flow times. After interpolating the two-dimensional 
axisymmetric solution onto the 3D grid, the shock cone 
adjusts slightly over a few flow times. By a time of 
AORa/Voo the mass accretion rate has settled to a con- 
stant value of 0.67 Mhl, as shown in Figured] This is 
co nsist ent with the average value of 0.56 Mhl reported 
by Ira], given the large variability seen in their models. 
In contrast, the mass accretion rate in our large accretor 
model varied by less than 0.02% once it settled down to 
the steady-state value. These small variations about the 
average could be attributed to numerical noise behind 



Figure 3. The pressure in a 2D plane a distance 0.8 R a down- 
stream of the accretor illustrates the highly axisymmetric nature 
of the downstream flow. The incoming flow is moving down the 
x-axis. A portion of the bow shock (on the Yin grid) is is also 
shown. The surface of the accreting star is shown with the leading 
face colored dark red to represent low mass flux and the back side 
colored yellow to represent high (inward) mass flux. (A color figure 
is available in the online journal.) 



The structure of this HLA flow is in reasonable agree- 
ment with time-independent, axisymmetric calculations 
IHunt] (|1971l ) . The local Mach number of the flow is shown 
in Figure 2] illustrating an upstream subsonic stagna- 
tion region, and a downstream, supersonic accretion flow. 
The region between the accretor surface and the stand- 
off bowshock is highly subsonic. The high pressure of 
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Figure 4. The Mach number in the post-shock flow for the large accretor, displayed on two orthogonal 
planes cutting through the flow axis. The surface of the accretor is colored by the local mass flux. Most of 
the accretion occurs on the downstream side of the accretor, at Mach numbers slightly greater than unity. 
(A color figure is available in the online journal.) 



this region results in a relatively low mass accretion rate 
through the upstream half of the accretor surface. In 
contrast, the inward accretion flow on the downstream 
side of the accretor reaches a Mach number of ~ 1.5. 
Most of the mass accretion occurs preferentially on the 
downstream side of the accretor; the local mass flux on 
the leading side is three times lower than that on the 
downstream side. At equilibrium, the bow shock has a 
stand-off distance of ~ 0.15 R a . While the opening angle 
of the bowshock is ill-defined in the vicinity of the ac- 
cretor, at a distance of several R a downstream the open- 
ing half-angle gradually decreases to a roughly constant 
~ 30°. This is noticeably larger than the theoretical 
prediction of sin _1 (l/M) » 20° for a Mach number of 
M = 3. (Petrich et. al. 1988 also noted a larger opening 
angle) . 

3.2. Small Accretor: R s /R a = 0.01 

Our second simulation reduced the accretor size to a 
value of R s /R a = 0.01, the same as in the T10 simulation 
of RA. At the time, RA's simulation lacked the duration 
to draw appropriate conclusions, although it displayed 
behavior that supported their hypothesized trend of more 
volatile behavior in M and J with decreasing accretor 
size. 

Our simulation ex tend s to 17 flow times, roughly three 
times longer than in IRAI and long enough to see the long 
term behavior unrelated to the initial conditions. The 
general structure of the accretion flow, illustrated in Fig- 
ure El remains similar to that seen in the larger accretor. 
The shock front is located at an upstream distance of 
~ 0.25 R a , significantly farther out than in the large- 
accretor simulation. The downstream opening half-angle 
is ~ 34°. Again the mass accretion occurs preferentially 
on the downstream side, with Mach numbers approach- 



ing 1.5. 

The most noticeable difference in this flow is the pres- 
ence of a breathing mode, in which the accretion wake os- 
cillates in the direction of the original flow while remain- 
ing axisymmetric. This oscillation begins immediately, 
and while it persists throughout the simulation, its am- 
plitude remains finite and the large-scale structure of the 
wake remains undisrupted. The stand-off distance of the 
bow shock varies quasi-periodically with an amplitude of 
about 10% of the steady-state value and an oscillation 
period slightly smaller than R a /Voo- The streaks in Fig- 
ure [5] betray the unsteady nature of the accretion flow 
onto the small accretor in contrast to the very steady 
flow onto the large accretor. 

The mass accretion rate, plotted in Figure El follows 
the same general evolution as for the large accretor. The 
accretion rate initially decreases until settling into a con- 
stant average after about ten flow times. In this case 
the late-time average is about 0.56 Mhl, somewhat lower 
than for the larger accretor. The main difference is an 
oscillation imposed on top of this long term behavior. 
These oscillations in the mass accretion rate have the 
same period as the oscillations in the shock stand-off 
distance, hence we attribute both to the axisymmetric 
breathing mode. The phase difference between mass ac- 
cretion and shock position is such that the brief minima 
in the mass accretion rate corresponds to when the lead- 
ing edge of the bowshock is falling back in toward the 
accretor. 

Although the flow is no longer steady, it remains highly 
axisymmetric. The accreted angular momentum is about 
an order of magnitude larger than in the large accretor 
simulations, but the values remain relatively small and 
there is no preferred orientation of the rotational flow. 
Even with the smaller radius of accretor surface, the in- 
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Figure 5. The Mach number in the post-shock flow for the small accretor, displayed on two orthogonal 
planes cutting through the flow axis. The surface of the accretor is colored by the local mass flux. (A color 
figure is available in the online journal.) 
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Figure 6. The mass accretion rate as a function of time for the 
model Rs/Ra = 0.01. 

stantaneous value of the specific angular momentum of 
accreting gas remains less than 5% of the specific angular 
momentum of a Keplerian orbit at the accretor surface. 

4. CONCLUSIONS 

We present high-resolution simulations of Hoyle- 
Lyttleton Accretion that demonstrate that the axisym- 
metry of HLA is remarkably stable. For both simulations 
we observe no rotational flow in the vicinity of the accre- 
tor, as indicated by the small amount of specific angular 
momentum in the accreted gas. For the large accretor 
case (R s — 0.05R a ) the mass accretion rate relaxes to 
a constant value of 0.67Mhx, with variations no greater 
than 0.02%. 

Decreasing the accretor size to R s = 0.01i? a introduces 
a quasi-periodic axisymmetric breathing mode. This 
mode appears almost immediately and reaches a quasi- 
stable oscillation within 10 R a /Voo. The mass accretion 
rate is modulated by the oscillation of the bow shock, 



but maintains an average of 0.56Mhl- The oscillations 
occur with a period of 0.86i? a /woo. This is decidedly 
larger than the acoustic period (the time for a sound wave 
to travel from the bowshock to the accretor and back), 
which we calculate to be 0.26^/^00. This breathing 
mode does not affect axisymmetry of the flow. 

Our results are dramatically different from those of 
RA, who used a similar numerical algorithm but a differ- 
ent grid geometry (nested cartesian grids versus overset 
spherical grids) with substantially lower spatial resolu- 
tion. Their simulations with R a /R s — 0.05 and 0.01 
both displayed quiescent periods broken by brief tran- 
sient disk-like flow. The larger accretor simulation in 
particular displayed an anti-correlation between the mass 
accretion rate and the specific angular momentum of the 
accreted gas: during the active periods when the accret- 
ing gas had a larger specific angular momentum there 
was a drop (by a factor of 2 or 3) in the mass accretion 
rate. This behavior was not as dramatic as the flip-flop 
flow observed in 2D HLA simulations. In contrast, our 
simulations do not display rotational flow and remain 
almost entirely axisymmetric throughout their runtimes. 

The simulations of RA included a random density per- 
turbation of 3% in their initial upstream flow as a means 
of breaking the initial symmetry of the flow. Our simula- 
tions do not include any initial noise, but some intrinsic 
noise is generated at the bow shock due to the numerical 
algorithm. The effects of this noise can be seen in Figure 
in which the asymmetry is visible near the downstream 
surface of the accretor. The addition of small density 
perturbations is likely to be insufficient to significantly 
affect the stability of the flow, and would likely have an 
effect simil ar to that of the nu merical noise. The sim- 
ulations of iRuff crt (1993, 11999( 1 suggest that large-scale 
density and velocity gradients (3% and 20% across the 
upstream accretion column) may be more likely to insti- 
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gate a non-stationary flow. These simulations displayed 
a relatively stable bow shock and no flip-flop behavior. 

The presence of the breathing mode and slight asym- 
metries near the accretor in the small-accretor case indi- 
cate that our results follow the general trend present in 
2D simulations: a smaller accretor produces a less stable 
flow pattern. This supports the notion that extending 
3D simulations to smaller accretor radii is critical for un- 
derstanding the true behavior of HLA in astrophysical 
contexts. The continued investigation of its stability un- 
der the influence of greater perturbations (including large 
scale gradients) and smaller accretor size will allow us to 
reassess what conditions, if any, produce an unstable flow 
and possible accretion of angular momentum. In addi- 
tion, introducing strong localized density perturbations 
analogous to those of a clumpy stellar wind may provide 
insight into the nature of supergiant fast x-ray transient 
outbursts while simultaneously addressing fundamental 
questions regarding stability. 
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